#creation d'une liste de SNPs à eliminer en fonction de hardy-weinberg, "minor allele frequency" et du pourcentage de missing
snp_selection <- function(inbfile=NULL,HW=NULL, miss=NULL, MAF=NULL,outputfile="SNPtoDel.txt"){
	if(!is.null(inbfile)){

		data <- NULL
		tmp_list <- NULL
		delHW <- NULL
		delmiss <- NULL
		delMAF <- NULL

		# creation d'un nombre aléatoire servant de prefixe
		randnb <- round(runif(10,min=0,max=1000))
		prefix <- randnb[1]

		cat("create SNPs list to be deleted : \n")

		if(!is.null(HW)){
			# creation d’un fichier de desequilibre de Hardy Weinberg (hwe.log, hwe.hh, hwe.hwe)
			system(paste("plink --noweb --bfile ",inbfile," --hardy --out ",prefix,"hwe",sep=""))
			hwe <- read.table(file=paste(prefix,"hwe.hwe",sep=""),h=T)
			
			# recherche du status des individus du fichier(cas, contrôle ou cas/contrôle)
			fam <- read.table(file=paste(inbfile,".fam",sep=""))
			fam2 <- fam[fam$V6 != -9,]
			study <- length(unique(fam2$V6)) # study = 1 (cases or controls), study = 2 (cases and controls) 

			if(study==1){
				if(unique(fam2$V6) == 1){hwe <- hwe[hwe$TEST=="UNAFF",];}
				if(unique(fam2$V6) == 2){hwe <- hwe[hwe$TEST=="AFF",];}

				del_hwe <- hwe[hwe$P < HW & !is.na(hwe$P),] 

				# creation de la liste de SNPs a supprimer
				delHW <- c(as.character(del_hwe$SNP))
			}
			else if(study == 2){
				aff <- hwe[hwe$TEST=="AFF",]
				unaff <- hwe[hwe$TEST=="UNAFF",]
				tmp <- merge(aff,unaff,by=c("CHR","SNP"))

				tmp_aff <- tmp[tmp$P.x < HW & !is.na(tmp$P.x),"SNP"] 
				tmp_unaff <- tmp[tmp$P.y < HW & !is.na(tmp$P.y),"SNP"]

				# creation de la liste de SNPs a supprimer
				delHW <- c(as.character(tmp_aff),as.character(tmp_unaff))
			}
			
			cat("add SNPs (p-value hwe < ",HW,")\n",sep="")
		}

		if(!is.null(miss)){
			# creation d'un fichier de difference de call rate (% de genotypage) entre cas et contrôles (dmiss.log, dmiss.hh, dmiss.missing)
			system(paste("plink --noweb --bfile ",inbfile," --test-missing --out ",prefix,"dmiss",sep=""))	
			dmiss <- read.table(file=paste(prefix,"dmiss.missing",sep=""),h=T)

			if(study==1){
				if(unique(fam2$V6) == 1){del_dmiss <- dmiss[dmiss$F_MISS_U > miss,"SNP"];}
				if(unique(fam2$V6) == 2){del_dmiss <- dmiss[dmiss$F_MISS_A > miss,"SNP"];}

				#  creation d'une liste de  SNPs avec un % de missing  > miss
				delmiss <- as.character(del_dmiss)
			}
			
			if(study==2){
				tmp_aff2 <- dmiss[dmiss$F_MISS_A > miss & !is.na(dmiss$F_MISS_A),"SNP"]
				tmp_unaff2 <- dmiss[dmiss$F_MISS_U > miss & !is.na(dmiss$F_MISS_U),"SNP"]

				# creation d'une liste de SNPs avec un % de missing  > miss
				delmiss <- c(as.character(tmp_aff2),as.character(tmp_unaff2))
			}

			cat("add SNPs (missing < ",miss,")\n",sep="")
		}

		if(!is.null(MAF)){
			# creation d'un fichier de frequence allelique (maf.log, maf.hh, maf.frq)
			system(paste("plink --noweb --bfile ",inbfile," --freq --out ",prefix,"maf",sep=""))
			maf <-  read.table(file=paste(prefix,"maf.frq",sep=""),h=T)

			# creation d'une liste de SNPs consideres comme "non frequent" (< MAF)
			delMAF <- maf[maf$MAF < MAF,"SNP"] 
			cat("add SNPs (minor allele frequency < ",MAF,")\n",sep="")
		}


		# concaténation des données et suppression des doublons
		delSNP <- c(delHW,delmiss,delMAF)
		final_list <- unique(delSNP)

		# enregistrement des donnees
		write.table(final_list,file=outputfile,quote=F, row.names=F)
		cat(paste("see results : ", outputfile,"\n",sep=""))

	}
	else { cat("missing inbfile parameter.\n")}
}


#####


######################### ANALYSE PAR PLAQUE #############################
###(creation des qqplots de comparaison de plaque les unes contre les autres et manhattan plot par plaque)###

# sert pour la fonction plate.analysis (renvoie le nom de la colonne qui correspond au minimum)
name <- function(data){
	ind <- which.min(data)
	val <- names(data)[ind]
}

# sert pour la fonction plate.analysis (calcul du lambda median)
lambda.median <- function(data){
	res <- (median(data,na.rm=T))/0.456
}

plate.analysis <- function(inbfile=NULL,plate_list=NULL,nb_by_graph=NULL,IC=NULL,outputfile="qqplot-plates"){

	if(is.null(inbfile) | is.null(plate_list) | is.null(nb_by_graph) | is.null(IC)){
		stop("ERROR : missing parameter. Please check : inbfile, plate_list, nb_by_graph and IC")
	}

	# creation d'un nombre aléatoire servant de prefixe
	randnb <- round(runif(10,min=0,max=1000))
	prefix <- randnb[1]

	system(paste("plink --noweb --bfile",inbfile,"--loop-assoc",plate_list,"--out",prefix,sep=" "))

	pat <- paste("^", prefix,".*.assoc$",sep="")

	# recuperation des fichiers contenant les informations de différence de fréquence allélique
	assocfiles <- dir(pattern=pat) # recuperation de tous les fichiers dont le nom contient le motif ".assoc" 
	p.names <-  NULL

	# chargement des donnees de plaque
	for ( i in 1:length(assocfiles)) {

		current_file <- read.table(file=assocfiles[i],h=T)

		# recuperation des colonnes : CHR, SNP et BP
		if(i == 1) { 
			pval <- current_file[c(1,2,3)]
			chisq <-current_file[c(1,2,3)]
		}

		tmp_pval <- current_file$P
		tmp_chisq <- current_file$CHISQ

		#  fusion des donnees
		pval <- cbind(pval,tmp_pval)
		chisq <- cbind(chisq,tmp_chisq)

		# recuperation du nom des plaques
		splitted_filename <- unlist(strsplit(assocfiles[i],".",fixed=TRUE)) # découpe le nom du fichier avec un "." comme delimiteur
		p.names <- c(as.character(p.names),as.character(splitted_filename[2]))
		write.table(p.names,file="plates-names.txt",quote=F,row.names=F,col.names=F)
	}

	#  on nomme les colonnes
	col_names <- c("CHR","SNP","BP",p.names)
	colnames(pval) <- col_names
	colnames(chisq) <- col_names
	
	pval_complete <- pval
	chisq_complete <- chisq

	# moyenne des chisq par plaque 
	chisq_tmp <- chisq[-c(1,2,3)]
	mean_chisq <- unlist(apply(X=chisq_tmp,na.rm=T , MARGIN=2, FUN=mean))

	# mediane des chisq par plaque
	lambda_median_chisq <- unlist(apply(X=chisq_tmp , MARGIN=2, FUN=lambda.median))

	# minimum des p-value par plaque
	pval_tmp <- pval[-c(1,2,3)]
	min_pval <- unlist(apply(X=pval_tmp,na.rm=T , MARGIN=2, FUN=min))
	
	# creation d'un fichier info_chisq_pval
	info_chisq_pval <- cbind(p.names,mean_chisq,lambda_median_chisq,min_pval)	

	# ajout de la p-value minimum 
	suppressWarnings(pval$P_min <- apply(X=pval_tmp, MARGIN=1, function(data) min(data,na.rm=TRUE)))

	# remplacement des NAs
	replaceNA <- apply(X=pval_tmp,MARGIN=1,function(data) ifelse(is.na(data),2,data))
	pval_tmp <- t(replaceNA)

	# ajout du nom de la plaque associee	
	P_name <- unlist(apply(X=pval_tmp, MARGIN=1, FUN=name))
	pval2 <- cbind(pval,P_name)

	# enregistrement des donnees
	datafile <- unlist(strsplit(inbfile,"/",fixed=TRUE))
	dataname <- datafile[length(datafile)]
	write.table(pval_complete,file=paste("p-values-complete-",dataname,".txt",sep=""),quote=F,row.names=F)
	write.table(chisq_complete,file=paste("chisq-complete-",dataname,".txt",sep=""),quote=F,row.names=F)

	write.table(pval2,file=paste("p-values-",dataname,".txt",sep=""), quote=F, row.names=F)
	cat(paste("see p-values-",dataname,".txt \n",sep=""))
	write.table(info_chisq_pval, file=paste("info_chisq_pval_",dataname,".txt",sep=""), quote=F, row.names=F,sep="\t")
	cat(paste("see info_chisq_pval_",dataname,".txt \n",sep=""))
	

	# creation des QQplots
	qqplot.p(dataname, p.names,nb_by_graph, IC, outputfile)
}


#####


# creation d'un QQplot par plaque - a faire pour les cas et pour les controles
qqplot.p<- function (dataname,p.names, nb_by_graph=4, IC=NULL, outputfile="qqplot-plates") #nb_by_graph : nb de plaque par graphique ; 
{
	if((nb_by_graph > length(p.names)) | (nb_by_graph==0)) {
		cat ("ERROR : incorrect nb_by_graph \n")
	}
	else{

		#  visualisation des p-values des plaques 
		nb_plate <- length(p.names)
		nb_plot <- ifelse(nb_plate%%nb_by_graph == 0, nb_plate%/%nb_by_graph, (nb_plate%/%nb_by_graph)+1)
		nb_page <- ifelse(nb_plot%%4 == 0, nb_plot%/%4, (nb_plot%/%4)+1)
		x<-1
		y<-nb_by_graph

		# chargement des donnees
		pval <- read.table(file=paste("p-values-",dataname,".txt",sep=""),h=T)
		p.names <- colnames(pval[c(4:(length(pval)-2))])
	}

	for (k in 1:nb_page){
			png(paste(outputfile,k,".png", sep=""), width=940, height=940, res=72)
			par(mfrow=c(2,2))		

			if(nb_page >1){
				if (k < nb_page){
					for (j in 1:4){ 
						plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
						nb_plot <- nb_plot-1
						x<-x+nb_by_graph
						y <- y+nb_by_graph
					}
				}
				else{
					if(y > length(p.names)){y<-length(p.names)}	
					for (i in 1:(nb_plot)){
						plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
						nb_plot <- nb_plot-1
						x <- x+nb_by_graph
						if((y+nb_by_graph) < length(p.names)){y <- y+nb_by_graph}
						else if(x <= length(p.names)){
							y <- length(p.names)
							plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
							nb_plot <- nb_plot-1
							break
							}
						}
					}
				}
			else{
				for (h in 1:nb_plot){ 
					plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
					x <- x+nb_by_graph
					if((y+nb_by_graph) < length(p.names)){y <- y+nb_by_graph}
					else if(x <= length(p.names)){
						y<- length(p.names)
						plot.QQs(data1=pval,pobs.name=p.names[x:y],IC=IC)
						break
					}
				}
			}
		dev.off()
	}
}


############################# FONCTION qqmultiple ##########################
plot.QQs <- function(	data1 = dat, 
			subset = NULL, #"ma_name == \"Hind\" & ma_type == \"affy\"",
			pobs.name = "pvalue", #c("pvalTest1", "pvalTest2")
			lim.X = NULL, ##c(xmin, xmax)
			lim.Y = NULL, ## c(ymin, ymax), 
			titre = NULL,
			couleur = T,	
			IC=NULL){

qqplot <- function (x, y, plot.it = TRUE, xlab = deparse(substitute(x)), 
    ylab = deparse(substitute(y)), add.to.plot = FALSE, ...) 
{
    sx <- sort(x)
    sy <- sort(y)
    lenx <- length(sx)
    leny <- length(sy)
    if (leny < lenx) 
        sx <- approx(1:lenx, sx, n = leny)$y
    if (leny > lenx) 
        sy <- approx(1:leny, sy, n = lenx)$y
    if (plot.it)
      if(add.to.plot)
        points(sx, sy, ...)
      else
        plot(sx, sy, xlab = xlab, ylab = ylab, ...)
    invisible(list(x = sx, y = sy))
}

confidence_intervals<- function (x, y)
{
    #x vector of pval
	#y confidence interval
	
	observed <- -(log10(x))
	N <- length(x) # number of p-values
	uniform <- -log10(1:N/(N+1))
	MAX <- max(c(observed, uniform))
	


	##Calcul d'interval de confiance
	
	q05 <- rep(0,N)
	q95 <- rep(0,N)
	
	for(i in 1:N)
	{
		#La k-ème statistique d'ordre d'un n-échantillon de lois uniformes U(0,1] suit la loi Beta(k,n-k+1) 
		q05[i] <- qbeta(1-y,i,N-i+1)
		q95[i] <- qbeta(y,i,N-i+1)
		
	}

	points(uniform,-log10(q95),col="black",type="l", lty=2, lwd=2)
	points(uniform,-log10(q05),col="black",type="l", lty=2, lwd=2)
	
}

print(dim(data1))
											
	if(!is.null(subset)) data1 <- data1[with(data1, eval(parse(text = subset))), ]
	
print(dim(data1))

numPval <- vector(length= length(pobs.name))

	#scale
	if(-log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T)) == -log10(0)) {
		min <- -log10(1e-200)
	}
	else{ 
		min <-  -log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T))
	}

#	if(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T))== -log10(1)) max <- 1
#		else max <-  -log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T))

#	if(is.null(lim.Y))  testY <- c(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T)), -log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T)))
	if(is.null(lim.Y))  testY <- c(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T)), min)
		else testY <- lim.Y

	if(is.null(lim.X)) testX <- testY
		else testX <- lim.X

for(i in 1:(length(pobs.name))){

 cat(paste("\n**********\ngraphe", i, ":", pobs.name[i], "\n**********\n"))

	pobs <- with(data1, eval(parse(text = pobs.name[i])))

	if(is.factor(pobs)) pobs <- as.numeric(levels(pobs))[pobs]

 cat(paste("nb p-values renseignees :", sum(!is.na(pobs))))

	numPval[i] <- sum(!is.na(pobs))

	#informations
	if(is.null(lim.X))	scaleX <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleX <- lim.X
	
	if(is.null(lim.Y))  scaleY <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleY <- lim.Y
	
cat("\nlimites X : "); print(scaleX)
cat("limites Y : "); print(scaleY)
	
	pobsmin <- min(pobs, na.rm = T)
	pobsmax <- max(pobs, na.rm = T)

cat(paste("pobsmin :", pobsmin))
cat(paste("\npobsmax :", pobsmax))	

cat(paste("\n-log10(pobsmin) :", -log10(pobsmin)))	
cat(paste("\n-log10(pobsmax) :", -log10(pobsmax), "\n"))

	qqplot( -log10((1:length(pobs)/(length(pobs) + 1))), #rang

		-log10(pobs), 
		main = paste("QQ-plot", subset, titre, sep = "\n"), 
		ylab="-log10(observed p-value)", 
		xlab="-log10(uniform p-value)", 
		xlim = testX,
  		ylim = testY , 
		pch = i,
		cex = 0.6,
		col = ifelse(couleur, i+1 ,1),
		add.to.plot = ifelse(i > 1, T, F)

)

	abline(0, 1, col="black", lwd=2)

}

# confidence interval
if (!is.null(IC))
{
		confidence_intervals( pobs, y=IC )  
}	

#legend
legend("bottomright", paste(pobs.name, " (N = ",  numPval, ")", sep = ""), pch = 1:length(pobs.name), col = 2:(length(pobs.name)+1))


}




############################# FONCTION qqmultiple222222 ##########################
plot.QQ2s <- function(	data1 = dat, 
			subset = NULL, #"ma_name == \"Hind\" & ma_type == \"affy\"",
			pobs.name = "pvalue", #c("pvalTest1", "pvalTest2")
			lim.X = NULL, ##c(xmin, xmax)
			lim.Y = NULL, ## c(ymin, ymax), 
			titre = NULL,
			couleur = T,	
			IC=NULL){

qqplot <- function (x, y, plot.it = TRUE, xlab = deparse(substitute(x)), 
    ylab = deparse(substitute(y)), add.to.plot = FALSE, ...) 
{
    sx <- sort(x)
    sy <- sort(y)
    lenx <- length(sx)
    leny <- length(sy)
    if (leny < lenx) 
        sx <- approx(1:lenx, sx, n = leny)$y
    if (leny > lenx) 
        sy <- approx(1:leny, sy, n = lenx)$y
    if (plot.it)
      if(add.to.plot)
        points(sx, sy, ...)
      else
        plot(sx, sy, xlab = xlab, ylab = ylab, ...)
    invisible(list(x = sx, y = sy))
}

confidence_intervals<- function (x, y)
{
    #x vector of pval
	#y confidence interval
	
	observed <- -(log10(x))
	N <- length(x) # number of p-values
	uniform <- -log10(1:N/(N+1))
	MAX <- max(c(observed, uniform))
	


	##Calcul d'interval de confiance
	
	q05 <- rep(0,N)
	q95 <- rep(0,N)
	
	for(i in 1:N)
	{
		#La k-ème statistique d'ordre d'un n-échantillon de lois uniformes U(0,1] suit la loi Beta(k,n-k+1) 
		q05[i] <- qbeta(1-y,i,N-i+1)
		q95[i] <- qbeta(y,i,N-i+1)
		
	}

	points(uniform,-log10(q95),col="black",type="l", lty=2, lwd=2)
	points(uniform,-log10(q05),col="black",type="l", lty=2, lwd=2)
	
}

											
	if(!is.null(subset)) data1 <- data1[with(data1, eval(parse(text = subset))), ]
	
print(dim(data1))

numPval <- vector(length= length(pobs.name))

	#scale
	if(-log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T)) == -log10(0)) min <- -log10(1e-200)
		else min <-  -log10(min(data1[,which(names(data1) %in% pobs.name)],na.rm=T))

	if(is.null(lim.Y))  testY <- c(-log10(max(data1[,which(names(data1) %in% pobs.name)],na.rm=T)), min)
		else testY <- lim.Y

	if(is.null(lim.X)) testX <- testY
		else testX <- lim.X

for(i in 1:(length(pobs.name))){

 cat(paste("\n**********\ngraphe", i, ":", pobs.name[i], "\n**********\n"))

	pobs <- with(data1, eval(parse(text=as.character(pobs.name[i]))))

	if(is.factor(pobs)) pobs <- as.numeric(levels(pobs))[pobs]

 cat(paste("nb p-values renseignees :", sum(!is.na(pobs))))

	numPval[i] <- sum(!is.na(pobs))

	#informations
	if(is.null(lim.X))	scaleX <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleX <- lim.X
	
	if(is.null(lim.Y))  scaleY <- c(-log10(max(pobs, na.rm = T)), -log10(min(pobs, na.rm = T)))
		else scaleY <- lim.Y
	
cat("\nlimites X : "); print(scaleX)
cat("limites Y : "); print(scaleY)
	
	pobsmin <- min(pobs, na.rm = T)
	pobsmax <- max(pobs, na.rm = T)

cat(paste("pobsmin :", pobsmin))
cat(paste("\npobsmax :", pobsmax))	

cat(paste("\n-log10(pobsmin) :", -log10(pobsmin)))	
cat(paste("\n-log10(pobsmax) :", -log10(pobsmax), "\n"))

	qqplot( -log10((1:length(pobs)/(length(pobs) + 1))), #rang

		-log10(pobs), 
		main = paste("QQ-plot", subset, titre, sep = "\n"), 
		ylab="-log10(observed p-value)", 
		xlab="-log10(uniform p-value)", 
		xlim = testX,
  		ylim = testY , 
		pch = i,
		cex = 0.6,
		col = ifelse(couleur, i+1 ,1),
		add.to.plot = ifelse(i > 1, T, F)

)

	abline(0, 1, col="black", lwd=2)

}

# confidence interval
if (!is.null(IC))
{
		confidence_intervals( pobs, y=IC )  
}	

#legend
legend("bottomright", paste(pobs.name, " (N = ",  numPval, ")", sep = ""), pch = 1:length(pobs.name), col = 2:(length(pobs.name)+1))


}
